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A criterion to identify the equilibration time in lipid bilayer simulations 

Rodolfo D. PorassoP* J.J. Lopez CascalesP 

With the aim of establishing a criterion for identifying when a lipid bilayer has reached 
steady state using the molecular dynamics simulation technique, lipid bilayers of different 
composition in their liquid crystalline phase were simulated in aqueous solution in pres- 
ence of CaCl2 as electrolyte, at different concentration levels. In this regard, we used two 
different lipid bilayer systems: one composed by 288 DPPC (DiPalmitoylPhosphatidyl- 
Choline) and another constituted by 288 DPPS (DiPalmitoylPhosphatidylSerine). In this 
sense, for both type of lipid bilayers, we have studied the temporal evolution of some lipids 
properties, such as the surface area per lipid, the deuterium order parameter, the lipid 
hydration and the lipid-calcium coordination. From their analysis, it became evident how 
each property has a different time to achieve equilibrium. The following order was found, 
from faster property to slower property: coordination of ions « deuterium order parameter 
> area per lipid m hydration. Consequently, when the hydration of lipids or the mean area 
per lipid are stable, we can ensure that the lipid membrane has reached the steady state. 



I. Introduction 

Over the last few decades, different computational 
techniques have emerged in different fields of sci- 
ence, some of them being extensively implemented 
and used by a great number of scientists around 
the globe. Among others, the Molecular Dynamics 
(MD) simulation is a very popular computational 
technique, which is widely used to obtain insight 
with atomic detail of steady and dynamic proper- 
ties in the fields of biology, physics and chemistry. 
In this regard, a critical aspect that must be iden- 
tified in all the MD simulations is related to the re- 
quired equilibration time to achieve a steady state. 
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This point is crucial in order to avoid simulation ar- 
tifacts that could lead to wrong conclusions. Cur- 
rently, with the increment of the computing power 
accessible to different investigation groups, much 
longer simulation trajectories are being carried out 
to obtain reliable information about the systems, 
with the purpose of approaching the time scale of 
the experimental phenomena. However, even when 
this fact is objectively desirable without further ob- 
jections, nowadays, much longer equilibration times 
are arbitrarily being required by certain reviewers 
during the revision process. From our viewpoint, 
this should be thoroughly revised due to the follow- 
ing two main reasons: first, because it results in a 
limiting factor in the use of this technique by other 
research groups which cannot access to very ex- 
pensive computing centers (assuming that authors 
provide enough evidence of the equilibration of the 
system). Second, to avoid wasting expensive com- 
puting time in the study of certain properties which 
do not require such long equilibration times, once 
the steady state of the system has been properly 
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identified. 

Phospholipid bilayers are of a high biological rele- 
vance, due to the fact that they play a crucial role in 
the control of the diffusion of small molecules, cell 
recognition, and signal transduction, among oth- 
ers. In our case, we have chosen the Phosphatidyl- 
Choline (PC) bilayer because it has been very well 
studied by MD simulations [TH7] and experimen- 
tally as well [8TH4], Furthermore, studies of the 
effects of different types of electrolytes on a PC bi- 
layer have also been studied, experimentally [T5"H2~4] 
and by simulation |25H33j . 

As mentioned above, the Molecular Dynamics 
(MD) simulations have emerged during the last 
decades as a powerful tool to obtain insight with 
atomic detail of the structure and dynamics of lipid 
bilayers [34 36 . Several MD simulations of mem- 
branes under the influence of different salt con- 
centrations have been carried out. One of the 
main obstacles related to these studies has been 
the time scale associated to the binding process of 
ions to the lipid bilayer. Considering the literature, 
a vast dispersion of equilibration times associated 
to the binding of ions to the membrane has been 
reported, where values ranging from 5 to 100 ns 
have been suggested for monovalent and divalent 
cations [25, 27-29, 32, 37, 38 . In this regard, we 
carried out four independent simulations of a lipid 
bilayer formed by 288 DPPC in aqueous solutions, 
for different concentrations of CaCl2 to provide an 
overview of their equilibration times. Among other 
properties, the surface area per lipid, the deuterium 
order parameters, lipid hydration and lipid-calcium 
coordination were studied. 

Finally, in order to generalize our results, a bi- 
layer formed by 288 DPPS in its liquid crystalline 
phase, in presence of CaCl2 at 0.25N, was simulated 
as well . 



II. Methodology 

Different molecular Dynamics (MD) simulations of 
lipid bilayer formed by 288 DPPC were carried out 
in aqueous solutions for different concentrations of 
CaCl2, from 0, up to 0.50 N. Furthermore, with 
the aim of generalizing our results, a bilayer of 288 
DPPS in presence of CaCl2 at 0.25 N was simulated 
as well. Note that the concentration of CaCl2 in 
terms of normality is defined as: 



Type of Lipid 


[CaCl 2 ] N 


Ca 2+ 


cr 


Water 


DPPC 











10068 
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0.06 
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10 


10053 


DPPC 


0.13 


12 


24 


10032 


DPPC 


0.25 


23 


46 


9999 


DPPC 


0.50 


46 


92 


9930 


DPPS 


0.25 


204 


120 


26932 



Table 1: The simulated bilayer systems. Note that 
the salt concentration is given in normal units. The 
numerals describe the number of molecules con- 
tained in the simulation box. 




DPPC 
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Figure 1: Structure and atom numbers for DPPC 
and DPPS used in this work. 
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equivalent weight = 
charge of the ions in the solution. 

In Table [TJ the number of molecules that consti- 
tute each system, applying Eq. (TTJ), is summarized. 

To build up the original system, a single DPPC 
lipid molecule, or DPPS lipid (Fig. [T]), was placed 
with its molecular axis perpendicular to the mem- 
brane surface [xy plane). Next, each DPPC, or 
DPPS, was randomly rotated and copied 144 times 
on each leaflet of the bilayer. Finally, the gaps ex- 
isting in the computational box (above and below 
the phospholipid bilayer) were filled using an equi- 
librated box containing 216 water molecules of the 
extended simple point charge (SPC/E) [32] water 
model. 

Thus, the starting point of the first system of 
Table [Q was formed by 288 DPPC in absence of 
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CaCl2- Once this first system was generated, the 
whole system was subjected to the steepest de- 
scent minimization process to remove any excess 
of strain associated with overlaps between neigh- 
boring atoms of the system. Thereby, the following 
DPPC systems in presence of CaCl2 were gener- 
ated as follows: to obtain a [CaCl2]=0.06 N, 15 
water molecules were randomly substituted by 5 
Ca 2+ and 10 CI - . An analogous procedure was ap- 
plied to the rest of the systems, where 36, 69 and 
138 water molecules were substituted by 12, 23 and 
46 Ca 2 + and 24, 46 and 92 CP, to obtain a [CaCl 2 ] 
concentration of 0.13 N, 0.25 N and 0.50 N, respec- 
tively. Finally, the DPPS bilayer was generated fol- 
lowing the same procedure described above for the 
DPPC, starting from a single DPPS molecule and 
once the lipid bilayer in presence of water passed 
the minimization process, 324 water molecules were 
substituted by 204 Ca 2+ and 120 Cl~ (note that 
144 of the 204 calcium ions were added to balance 
the negative charge associated with the DPPS). 

GROMACS 3.3.3 package [30HT] was used in the 
simulations, and the properties showed in this work 
were obtained using our own code. The force field 
proposed by Egberts et al. [2^ was used for the 
lipids, and a time step of 2 fs was used as integra- 
tion time in all the simulations. A cut-off of 1.0 
nm was used for calculating the Lennard- Jones in- 
teractions. The electrostatic interaction was evalu- 
ated using the particle mesh Ewald method [4"2"ll4"3"] . 
The real space interaction was evaluated using a 
0.9 nm cut-off, and the reciprocal space interaction 
using a 0.12 nm grid with a fourth-order spline in- 
terpolation. A semi-isotropic coupling pressure was 
used for the coupling pressure bath, with a refer- 
ence pressure of 1 atm which allowed the fluctu- 
ation of each axis of the computer box indepen- 
dently. For the DPPC bilayer, each component of 
the system (i.e., lipids, ions and water) was coupled 
to an external temperature coupling bath at 330 K, 
which is well above the transition temperature of 
314 K [HIES]. For DPPS bilayer, each component 
of the system was coupled to an external temper- 
ature coupling bath at 350 K, which is above the 
transition temperature (3S1I37]. All the MD sim- 
ulations were carried out using periodic boundary 
conditions. The total trajectory length of each sim- 
ulated system was of 80 ns of MD simulation, where 
the coordinates of the system were recorded every 
5 ps for their appropriate analysis. 



Finally, in order to study the effect of the temper- 
ature, only the case corresponding to 0.25 N CaCl2 
was investigated at two additional temperatures, 
340 K and 350 K. 



III. Results and discussion 

i. Effect of the CaCl2 concentration 

a. Surface area per lipid 

Surface area per lipid (A) is a property of lipid 
bilayers which has been accurately measured from 
experiments [48] . The calculation of mean area per 
lipid can be determined from the MD simulation 
as: 



(A) 



x ■ y 

N 



(2) 



where x and y represent the box sizes in the direc- 
tion x and y (perpendicular to the membrane sur- 
face) over the simulation, and N is the number of 
lipids contained in one leaflet, in our case N = 144. 

Focusing on the study of the time evolution of 
the area per lipid, Figure [2] depicts the running 
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Figure 2: Running area per lipid at T = 330 K 
in presence of [CaCl 2 ] at (A) 0.06 N, (B) 0.13 N, 
(C) 0.25 N, (D) 0.50 N and (E) 0.25 N (in this 
case, T = 350 K). Solid lines represent the mean 
area obtained from the last 70 ns of the simulated 
trajectories (see text for further explanation). The 
type of lipid is indicated in the legends. 



040005-3 



Papers in Physics, vol. 4, art. 040005 (2012) / R D Porasso et al. 



surface area per lipid for different concentrations 
of CaCi2 and type of lipid. In general, for the 5 
bilayers formed by DPPC or DPPS, the area per 
lipid achieved a steady state after 10 ns of simula- 
tion, being this equilibration time almost indepen- 
dent of the concentration of CaCl2 and type of lipid 
which composed the membrane. 

In absence of salt, an average area per lipid of 
(A) = 0.663 ± 0.008 nm 2 was calculated from the 
last 70 ns of the simulated trajectory, discarding 
the first 10 ns corresponding to the equilibration 
time. This value agrees with experimental data, 
where values in a range from 0.55 to 0.72 nm 2 have 
been measured [TOIIHl ligHST]. Table H shows the 
mean surface area per lipid (again, after discarding 
the equilibration time of 10 ns) with their corre- 
sponding error bar. 

From the simulation results, a shrinking in the 
surface area per lipid with the increment of the 
ionic strength of the solution is observed. This 
shrinking is expected and attributed to the com- 
plexation of lipid molecules by calcium, such as it 
has been pointed out in previous studies [28,29,52 . 

b. Deuterium order parameter 

The deuterium order parameter, Scd, is measured 
from 2 H-NMR experiments. This parameter pro- 
vides relevant information related to the disorder of 
the hydrocarbon region in the interior of the lipid 
bilayers by measuring the orientation of the hydro- 
gen dipole of the methylene groups with respect 
to the perpendicular axis to the lipid bilayer. Due 
to the fact that hydrogens of the lipid methylene 
groups (CH2) have not been taken into account (in 
an explicit way) in our simulations, the order pa- 
rameter —Scd on the i + 1 methylene group was 
defined as the normal unitary vector to the vector 
defined from the i to the i + 2 CH2 group and con- 
tained in the plane formed by the methylene groups 
i, i ± 1 and i + 2. Thus, the deuterium order pa- 
rameter —Scd on the i — th of the CH2 group can 
be estimated by Molecular Dynamics simulations 
as follows: 

- Scd = ~ (3 cos 2 (6) - 1) (3) 

where 9 is the angle formed between the unitary 
vector defined above and the z axis. The expres- 
sion in brackets (...) denotes an average over all 
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[CaCl 2 ] N 


(A) (nm 2 ) 


Hydration 


Lipid 






Number 


DPPC 





0.663 ±0.008 


1.758 ±0.009 


DPPC 


0.06 


0.658 ±0.008 


1.740 ±0.009 


DPPC 


0.13 


0.651 ±0.007 


1.719 ±0.010 


DPPC 


0.25 


0.641 ±0.009 


1.680 ±0.015 


DPPC 


0.50 


0.628 ±0.010 


1.610 ±0.015 


DPPS 


0.25 


0.522 ±0.007 


2.552 ±0.010 



Table 2: Area per lipid and lipid hydration number 
as a function of salt concentration (see text for fur- 
ther explanation) . Note that the salt concentration 
is given in normal units. Error bars were calculated 
for each system separately from subtrajectories of 
10 ns length. Simulation temperature = 330 K. 

the lipids and time. Hence, note that the —Scd 
can adopt any value between -0.5 (corresponding 
to a parallel orientation to the lipid/water inter- 
face) and 1 (oriented along the normal axis to the 
lipid bilayer). 

Figure |3] shows the running —Scd for different 
carbons of the DPPC and DPPS tails and salt con- 
centrations. Only the carbons which correspond to 
the initial (hydrocarbons 2 and 6), the middle (hy- 
drocarbon 10) and final (hydrocarbons 13 and 15) 
methylene groups of the lipid tails were depicted 
in this figure. Each point of the figure represents 
the average values of —Scd over 5 ns of subtrajec- 
tory length, and the lines represent the mean val- 
ues calculated from the last 70 ns of the trajectories 
simulated. From this figure, it is observed how in 
all the cases, the required equilibration time is less 
than 10 ns of simulation time, independently of the 
salt concentration and the type of lipid. Finally, it 
is noted that Figure [3] exhibits an increase in the 
deuterium order parameters with the salt concen- 
tration, consistent with the shrinking of the area 
per lipid described above. 

c. Lipid hydration 

To analyze the lipid hydration, the radial distri- 
bution function g(r) of water around one of the 
oxygens of the phosphate group (atom number 10 
in Fig. Q] for DPPC and DPPS) was calculated. 
The radial distribution function g (r) is defined as 
follows: 
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Figure 3: Running deuterium order parameter, 
-S C D, in presence of [CaCl 2 ] at (A) 0.06 N, (B) 
0.13 N, (C) 0.25 N, (D) 0.50 N and (E) 0.25 N. 
DPPC simulations were performed at 330 K and 
DPPS simulation temperature was 350 K. Solid 
lines represent mean values —Scd obtained from 
the last 70 ns of the simulated trajectories. The 
type of lipid is indicated in the legends. Symbols: 
o hydrocarbon 2; o hydrocarbon 6; <l hydrocarbon 
10; + hydrocarbon 13 and x hydrocarbon 15. Note 
that the error bars have the same size as symbol. 



g( r ) = 



N{r) 



(4) 



Aitr 2 pSr 

where N (r) is the number of atoms in a spherical 
shell at distance r and thickness 6r from a reference 
atom, p is the density number taken as the ratio of 
atoms to the volume of the total computing box. 

From numerical integration of the first peak of 
the radial distribution function, the hydration num- 
bers can be estimated for different atoms of the 
DPPC or DPPS. Figure S] depicts the hydration 
number of phosphate oxygen (atom 10 in Fig. [1] 
for DPPC and DPPS) in presence of CaCl2, where 
each point represents the average of 5 ns subtrajec- 
tory length. These results show how this property 
reached a steady state for the cases (A), (B) and 
(E), after 10 ns of simulation. However, for the 
cases (C) and (D), 5 ns of extra simulation trajec- 
tory were required to reach a steady state. Tabled] 
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time (ns) 



Figure 4: Hydration number of the phosphate oxy- 
gen (atom 10 in Fig. [1]) along the simulated tra- 
jectories in presence of [CaCy at (A) 0.06 N, (B) 
0.13 N, (C) 0.25 N, (D) 0.50 N (for DPPC T = 
330 K) and (E) 0.25 N. In this case, T = 350 K. 
Solid lines represent the mean value of the hydra- 
tion number calculated from the last 70 ns of the 
simulated trajectories. The type of lipid is indi- 
cated in the legends. Note that the error bars have 
the same size as symbol. 



shows the hydration numbers for the last 70 ns of 
the trajectory length, corresponding to four concen- 
trations of CaCi2 and both types of lipids, DPPC 
and DPPS. In this regard, from Fig. |4j the signif- 
icant lipid dehydration with the increment of the 
ionic strength of the solution is evident, in good 
accordance with previous results [5 2) . 

d. Phospholipid- calcium coordination 

Some authors have reported how the lipid coordi- 
nation by divalent cations widely varies . Thus, on 
the one hand, some authors [25] have reported that 
this is a very slow process, which requires about 85 
ns of simulation time, but, on the other hand, other 
authors [53] have suggested that this process results 
much more rapid, taking less than 1 ns. In this 
sense, the coordination of DPPC-Ca 2+ was studied 
by monitoring the oxygen-calcium coordination of 
the carbonyl oxygens (atoms 16 and 35 in Fig. [1]) 
and phosphate oxygens (atoms 9 and 10 of DPPC 
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Figure 5: Left column represents the number of 
Ca 2+ coordinated to lipids in presence of [CaCy 
at (A) 0.06 N, (B) 0.13 N, (C) 0.25 N and (D) 0.50 
N, for T = 330 K. Right column shows the quantity 
of calcium ions coordinated to lipids expressed in 
percentage, along the simulated trajectory. 

in Fig. [!}, as a function of time. The left column 
in Fig. [3] represents the oxygen coordination num- 
ber, while the right one depicts the percentage of 
calcium ions involved in the coordination process 
with respect to the total number of calcium ions 
present in the aqueous solution. Figure [5] shows 
how the DPPC coordination by calcium is a quick 
process, taking less than 5 ns of simulation time 
to achieve a steady state. The kinetic of this pro- 
cess appears to be related to the ratio between cal- 
cium/lipid. Thus, after the first 5 ns of simulation 
time, the Ca-lipid coordination presents some fluc- 
tuation along the rest of the simulated trajectory. 
However, in Fig. [5] (A) and (B) (for the cases of 
lower concentration), it is observed how the per- 
centage of coordination fluctuates between a 60% 
and a 100%. We consider that this broad fluctu- 
ation is related to the limited sample size of our 
simulations that introduces a certain noise in our 
results. 

ii. Effect of temperature 

This section focuses on the study of the role played 
by temperature on the equilibration process. In 
this regard, only the system corresponding to a 



T (K) {A) (nm 2 ) Hydration Number 

330 0.642 ±0.009 1.680 ± 0.010 

340 0.650 ±0.007 1.683 ± 0.020 

350 0.666 ±0.008 1.689 ± 0.015 



Table 3: Area per lipid and the lipid hydra- 
tion number as a function of temperature. Error 
bars were calculated from subtrajectories of 10 ns 
length. DPPC bilayer in presence of [CaCy = 0.25 
N. 



concentration of 0.25 N in CaCk) was studied, for 
a range of temperatures from 330 K to 350 K (all 
of them above the transition temperature of 314 K 
HUES] for the DPPC). 

Figure |5] shows the running area along the trajec- 
tory. In this case, it was noticed how the systems 
achieve a steady state after a trajectory length of 
roughly 10 ns, where Table [3] shows the mean area 
per lipid calculated from the last 70 ns of simu- 
lation time. Figure [7] shows the deuterium order 
parameter of the methylene groups along the lipid 
tails, calculated from Eq. Figure on the 

one hand, clearly shows that for the three temper- 
atures the systems have reached the steady state 
before the first 10 ns of simulation. On the other 
hand, it shows an increase in the disorder of the 
lipid tails with temperature, which is closely re- 
lated with the increase of the area per lipid, such 
as it was pointed out above. Figure |8] depicts the 
results of the hydration numbers of DPPC for the 
three temperatures studied, where the equilibrated 
state was achieved after 10 ns of simulation time. 
Table [3] provides the calculated hydration numbers 
in the equilibrium, showing how the lipid hydra- 
tion remained invariable with the rising of the tem- 
perature. Concerning the lipid-calcium coordina- 
tion, Fig. IHl represents the lipid-calcium coordina- 
tion number, and the right column represents the 
calcium that participates in the coordination ex- 
pressed in percentage respect the total of calcium 
ions in solution. From simulation, it becomes ev- 
ident how calcium ions required less than 5 ns to 
achieve an equilibrated state for the three temper- 
atures studied. In summary, for all the properties 
studied in this section, a slight decrease in the equi- 
libration time with the increasing temperature was 
appreciated. 
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Figure 6: Running area per lipid for [CaCy = 0.25 
N at different temperatures, (A) T = 330 K, (B) T 
= 340 K and (C) T = 350 K. Solid lines represent 
the mean values obtained from the last 70 ns of 
simulation. 
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Figure 7: Deuterium Order Parameter, —Scd, 
along the simulated trajectory for a concentration 
of [CaCy = 0.25 N, for the following temperatures: 
(A) T = 330 K, (B) T = 340 K and (C) T = 350 
K. Solid lines represent the average order parame- 
ter for the last 70 ns of simulation. Note that the 
error bars have the same size as symbol. 



IV. Conclusions 

The present work deals with the simulation time re- 
quired to achieve the steady state for a lipid bilayer 
system in presence of CaCy In this regard, wc 
studied two different systems: one with DPPC and 
another one with DPPS bilayer; both systems in 



Figure 8: Hydration number of phosphate oxygen 
(atom 10 in Fig. [TJ along the simulated trajectories 
for a [CaCy = 0.25 N at different temperatures: 
(A) T = 330 K, (B) T = 340 K and (C) T = 350 K. 
Solid lines represent the average hydration number 
for the last 70 ns of simulation. Note that the error 
bars have the same size as symbol. 

presence of CaCl2 (at different level concentration) . 
The salt free case was also studied, as control. The 
analysis of various lipid properties studied here in- 
dicates that some properties reach the steady state 
more quickly than others. In this sense, we found 
that the area per lipid and the hydration number 
are slower than the deuterium order parameter and 
the coordination of cations. Consequently, to en- 
sure that a system composed by a lipid bilayer has 
reached a steady state, the criterion that we pro- 
pose is to show that the area per lipid or the hy- 
dration number have reached the equilibrium. 

From our results, two important aspects should 
be remarked: 

1. The equilibration time is strongly depen- 
dent on the starting conformation of the sys- 
tem. Wrong starting conformations will re- 
quire much longer equilibration times, even of 
one order of magnitude higher than the re- 
quested from a more refined starting confor- 
mation. 

2. Temperature is a critical parameter for reduc- 
ing the equilibration time in our simulations, 
due to the fact that higher temperatures in- 
crease the kinetic processes, i.e., the sampling 
of the configurational space of the system. 
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Figure 9: Left column represents the number of cal- 
cium ions involved in the lipid coordination along 
time for a concentration of [CaCl2]= 0.25 N at dif- 
ferent temperatures: (A) T = 330 K, (B) T = 340 
K and (C) T = 350 K. The right column shows the 
same information expressed as a percentage of the 
total number of calcium ions in solution. 



Acknowledgements - Authors wish to thank the 
assistance of the Computing Center of the Universi- 
dad Politecnica de Cartagena (SAIT), Spain. RDP 
is member of 'Carrera del Investigador', CON- 
ICET, Argentine. 



[1] 



O Berger, O Edholm, F Jahnig, Molecular dy- 
namics simulations of a fluid bilayer of di- 
palmitoylphosphatidylcholine at full hydration, 
constant pressure, and constant temperature, 
Biophys. J. 72, 2002 (1997). 



[2] E Egberts, S J Marrink, H J C Berendsen, 
Molecular dynamics simulation of a phospho- 
lipid membrane, Eur. Biophys. J. 22, 423 
(1994). 

[3] U Essmann, L Perera, M L Bcrkowitz, The 
origin of the hydration interaction of lipid bi- 
layers from MD simulation of dipalmitoylphos- 
phatidylcholine membranes in gel and crys- 
talline phases, Langmuir 11, 4519 (1995). 



[4] S E Feller, Y Zhang, R W Pastor, R B Brooks, 
Constant pressure molecular dynamics simula- 
tions: The Langevin piston method, J. Chem. 
Phys. 103, 4613 (1995). 

[5] W Shinoda, T Fukada, S Okazaki, I Okada, 
Molecular dynamics simulation of the dipalmi- 
toylphosphatidylcholine (DPPC) lipid bilayer 
in the fluid phase using the Nosr-Parrinello- 
Rahman NPT ensemble, Chem. Phys. Lett. 
232, 308 (1995). 

[6] D P Tieleman, H J C Berendsen, Molecular dy- 
namics simulations of a fully hydrated dipalmi- 
toylphosphatidylcholine bilayer with different 
macroscopic boundary conditions and param- 
eters, J. Chem. Phys. 105, 4871 (1996). 

[7] K Tu, D J Tobias, M L Klein, Constant 
pressure and temperature molecular dynamics 
simulation of a fully hydrated liquid crystal 
phase dipalmitoylphosphatidylcholine bilayer, 
Biophys. J. 69, 2558 (1995). 

[8] M F Brown, Theory of spin-lattice relaxation 
in lipid bilayers and biological membranes. 
Dipolar relaxation, J. Chem. Phys. 80, 2808 
(1984). 

[9] M F Brown, Theory of spin-lattice relaxation 
in lipid bilayers and biological membranes. 
2 H and N quadrupolar relaxation, J. Phys. 
Chem. 77, 1576 (1982). 

[10] J F Nagle, R Zang, S Tristam-Nagle, W S Sun, 
H I Petrache, R M Suter, X-ray structure de- 
termination of fully hydrated L. phase dipalmi- 
toylphosphatidylcholine bilayers, Biophys. J. 
70, 1419 (1996). 

[11] R P Rand, V A Parsegian, Hydration forces be- 
tween phospholipid bilayers, Biochim. Biophys. 
Acta 988, 351 (1989). 

[12] J Seelig, Deuterium magnetic resonance: The- 
ory and application to lipid membranes, Q. 
Rev. Biophys. 10, 353 (1977). 

[13] J Seelig, A Seelig, Lipid conformation in model 
membranes and biological systems, Q. Rev. 
Biophys. 13, 19 (1980). 



040005-8 



Papers in Physics, vol. 4, art. 040005 (2012) / R D Porasso et al. 



[14] W J Sun, R M Sutcr, M A Knewtson, C R 
Worthington, S Tristram-Nagle, R Zhang, J F 
Nagle, Order and disorder in fully hydrated un- 
oriented bilayers of gel phase dipalmitoylphos- 
phatidylcholine, Phys. Rev. E. 49, 4665 (1994). 

[15] H Akutsu, J Seelig, Interaction of metal ions 
with phosphatidylcholine bilayer membranes, 
Biochemistry 20, 7366 (1981). 

[16] M G Gancsan, D L Schwinke, N Weiner, Effect 
of Ca 2+ on thermotropic properties of satu- 
rated phosphatidylcholine liposomes, Biochim. 
Biophys. Acta 686, 245 (1982). 

[17] L Herbette, C A Napolitano, R V McDaniel, 
Direct determination of the calcium profile 
structure for dipalmitoyllecithin multilayers 
using neutron diffraction, Biophys. J. 46, 677 
(1984). 

[18] D Huster, K Arnold, K Gawrisch, Strength of 
Ca 2+ binding to retinal lipid membrane: Con- 
sequences for lipid organization, Biophys. J. 
78, 3011 (2000). 

[19] Y Inoko, T Yamaguchi, K Furuya, T Mitsui, 
Effects of cations on dipalmitoyl phosphatidyl- 
choline/ cholesterol/ water systems, Biochim. 
Biophys. Acta 413, 24 (1975). 

[20] R Lehrmann, J J Seelig, Adsorption of Ca 2+ 
and La 3+ to bilayer membranes: Measurement 
of the adsorption enthalpy and binding con- 
stant with titration calorimetry, Biochim. Bio- 
phys. Acta 1189, 89 (1994). 

[21] L J Lis, W T Lis, V A Parsegian, R P Rand, 
Adsorption of divalent cations to a variety of 
phosphatidylcholine bilayers, Biochemistry 20, 
1771 (1981). 

[22] L J Lis, V A Parsegian, R P Rand, Binding 
of divalent cations to dipalmitoylphosphatidyl- 
choline bilayers and its effect on bilayer inter- 
action, Biochemistry 20, 1761 (1981). 

[23] T Shibata, Pulse NMR study of the interaction 
of calcium ion with dipalmitoylphosphatidyl- 
choline lamellae, Chem. Phys. Lipids. 53, 47 
(1990). 



[24] S A Tatulian, V I Gordeliy, A E Sokolova, A G 
Syrykh, A neutron diffraction study of the in- 
fluence of ions on phospholipid membrane in- 
teractions, Biochim. Biophys. Acta 1070, 143 
(1991). 

[25] R A Bockmann, H Grubmiillcr, Multistep 
binding of divalent cations to phospholipid bi- 
layers: A molecular dynamics study, Ange- 
wandte Chemie 43, 1021 (2004). 

[26] J Faraudo, A Travesset, Phosphatidic acid do- 
mains in membranes: Effect of divalent coun- 
terions, Biophys. J. 92, 2806 (2007). 

[27] A A Gurtovenko, Asymmetry of lipid bilay- 
ers induced by monovalent salt: Atomistic 
molecular- dynamics study, J. Chem. Phys. 
122, 244902 (2005). 

[28] P Mukhopadhyay, L Monticelli, D P Tieleman, 
Molecular dynamics simulation of a palmitoyl- 
oleoyl phosphatidyls erine bilayer with Na + 
counterions and NaCl, Biophys. J. 86, 1601 
(2004). 

[29] S A Pandit, D Bostick, M L Bcrkowitz, Molec- 
ular dynamics simulation of a dipalmitoylphos- 
phatidylcholine bilayer with NaCl, Biophys. J. 
84, 3743 (2003). 

[30] U R Pedersen, C Laidy, P Westh, G H Pe- 
ters, The effect of calcium on the properties 
of charged phospholipid bilayers, Biochim. Bio- 
phys. Acta 1758, 573 (2006). 

[31] J N Sachs, H Nanda, H I Petrache, T B Woolf, 
Changes in phosphatidylcholine headgroup tilt 
and water order induced by monovalent salts: 
Molecular dynamics simulations, Biophys. J. 
86, 3772 (2004). 

[32] K Shinoda, W Shinoda, M Mikami, Molec- 
ular dynamics simulation of an archeal lipid 
bilayer whit sodium chloride, Phys. Chem. 
Chem. Phys. 9, 643 (2007). 

[33] N L Yamada, H Seto, T Takeda, M 
Nagao, Y Kawabata, K Inouc, SAXS, 
SANS and NSE studies on unbound state 
in DPPC/water/CaCh system, J. Phys. Soc. 
Jpn. 74, 2853 (2005). 



040005-9 



Papers in Physics, vol. 4, art. 040005 (2012) / R D Porasso et al. 



[34] D Frenkel, B Smit, Understanding molecu- 
lar simulations, Academic Press, New York 
(2002). 

[35] J J Lopez Cascales, J Garcia de la Torre, S J 
Marrink, H J C Berendsen, Molecular dynam- 
ics simulation of a charged biological mem- 
brane, J. Chcm. Phys. 104, 2713 (1996). 

[36] W F van Gunstcrcn, H J C Berendsen, 
Computer simulations of molecular dynamics: 
Methodology, applications and perspectives in 
chemistry, Angew. Chcm Int. Ed. Engl. 29, 
992 (1990). 

[37] R A Bockmann, A Hac, T Heimburg, H 
Grubmiillcr, Effect of sodium chloride on a 
lipid bilayer, Biophys. J. 85, 1647 (2003). 

[38] A A Gurtovenko, I Vattulainen, Effect of NaCl 
and KCl on phosphatidylcholine and phos- 
phatidylethanolamine lipid membranes: In- 
sight from atomic-scale simulations for un- 
derstanding salt-induced effects in the plasma 
membrane, J. Phys. Chem. B. 112, 1953 
(2008). 

[39] H J C Berendsen, J R Grigera, T P Straatsma, 
The missing term in effective pair potentials, 
J. Phys. Chem. 91, 6269 (1987). 

[40] H J C Berendsen, D van der Spoel, R van 
Drunen, A message-passing parallel molecu- 
lar dynamics implementation, Comp. Phys. 
Comm. 91, 43 (1995). 

[41] E Lindahl, B Hess, D van der Spoel, GRO- 
MACS 3.0: A package for molecular simula- 
tion and trajectory analysis, J. Mol. Mod. 7, 
306 (2001). 

[42] T Darden, D York, L Pedersen, Particle mesh 
Ewald: An N.log(N) method for Ewald sums 
in large systems, J. Chcm. Phys. 98, 10089 
(1993). 

[43] U Essmann, L Perea, M L Berkowitz, T Dar- 
den, H Lee, L G Pedersen, A smooth particle 
mesh Ewald method, J. Chem. Phys. 103, 8577 
(1995). 



[44] L R De Young, K A Dill, Solute partition- 
ing into lipid bilayer-membranes, Biochem- 
istry 27, 5281 (1988). 

[45] A Seelig, J Seelig, The dynamic structure of 
fatty acyl chains in a phospholipid bilayer mea- 
sured by deuterium magnetic resonance, Bio- 
chemistry 13, 4839 (1974). 

[46] G Cevc, A Watts, D Marsh, Titration of 
the phase transition of phosphatidilserine bi- 
layer membranes. Effect of pH, surface electro- 
statics, ion binding and head-group hydration, 
Biochemistry 20, 4955 (1981). 

[47] H Hauser, F Paltauf, G G Shipley, Structure 
and thermotropic behavior of phosphatidyls er- 
ine bilayer membranes, Biochemistry 21, 1061 
(1982). 

[48] J F Nagle, S Tristam-Nagle, Structure of lipid 
bilayers, Biochim. Biophy. Acta 1469, 159 
(2000). 

[49] B A Lewis, D M Engelman, Lipid bilayer thick- 
ness varies linearly with acyl chain length in 
fluid phosphatidylcholine vesicles, J. Mol. Biol. 
166, 211 (1983). 

[50] R J Pace, S I Cham, Molecular motions in 
lipid bilayer. I. Statistical mechanical model of 
acyl chains motion, J. Chem. Phys. 76, 4217 
(1982). 

[51] R L Thurmond, S W Dodd, M F Brown, 
Molecular areas of phospholipids as deter- 
mined by 2H NMR spectroscopy, Biphys. J. 59, 
108 (1991). 

[52] R D Porasso, J J Lopez Cascales, Study of the 
effect of Na + and Ca 2+ ion concentration on 
the structure of an asymmetric DPPC'/DPPS 
+ DPPS lipid bilayer by molecular dynamics 
simulation, Coll. and Surf. B. Bioint. 73, 42 
(2009). 



040005-10 



